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MHD instabilities in accretion mounds - 1: 2D axisymmetric 
simulations 
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ABSTRACT 

We have performed stability analysis of axisymmetric accretion mounds on neutron stars in 
High Mass X-ray Binaries (HMXB) by 2-D MHD simulations with the PLUTO MHD code. 
We find that the mounds are stable with respect to interchange instabilities, but addition of 
excess mass destabilizes the equilibria. Our simulations confirm that accretion mounds are 
unstable with respect to MHD instabilities beyond a threshold mass. We investigate both filled 
and hollow mounds and for the latter also compute the expected profile of cyclotron resonance 
scattering features (CRSF). In comparison to the CRSF from filled mounds reported in our 
earlier work, hollow mounds display wider and more complex line profiles. 

Key words: accretion — magnetic fields — (stars:) binaries: general — X-rays: binaries — 
line: formation — radiation mechanisms: non-thermal 



1 INTRODUCTION 

Neutron stars in accreting X-ray pulsars ac crete matter from the 
comp anion star either from stellar winds (iDavidson & Ostrikerl 
1 19731) or through disc accretion by Roche lobe overflow 
dGhosh et all 1 19771 : iKoldoba et alJ [2002 : iRomanova et alJ l2003h . 
They can be broadly classified into two classes: 1) high mass X-ray 
binaries (HMXB) with companion stars of masses several times 
the solar mass and neutron stars with high surface magnetic field 
- 10 12 G and 2) low mass X-ray binaries (LMXB) with compan- 
ion stars of masses less than a solar mass and neutron star mag- 
netic fields several orders lower in mag nitude ~ 10 7 G — 10 9 G (see 
Bhattacharya & van den Heuvel (1991) for a review). In this paper 
we consider the effect of accretion on the evolution of surface mag- 
netic field of HMXB sources by the formation of accretion mounds. 

The accreted matter in HMXB passes through a shock, 
gradually settling down on the polar cap to form an accretion 
mound. X-ray emission from such mounds sho w characteristic cy- 
clotron resonance scatte r ing fe a tures (CRSF ) (Hardin g & Preecd 
19871: lArava & Harding! Il999l: lArava-Gochez & Hardingl l2000l: 
Becker & Wolfd l2007» . The CRSF depends on the magnetic field 
of the local emitting region, and hence serve as a tool to understand 
the structure of accretion columns. CRSF often show complex line 
features and characteristic variat ions with rotation phase and the 
luminosity of the neutron star dCoburn et al.l 120021: iHeindl et all 
|2004 iMihara et al.ll2007l : iLutovinov & Tsvgankovll2008l) . ExplairT 
ing such features require appropriate modelling of the structure of 
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the accretion column and the effect of accretion induced field dis- 
tortion from the accretion mound. 

Also, several authors propose that diamagnetic screen- 
ing of the magnetic field c an lower the apparent dipole 



moment of the neutron star dRomanil Il990l: Cumming et al. l 
200ll: iMelatos & Phinnev) l200ll : IChoudhuri & Konaj 120021 : 



Kona r & ChoudhuriT 120041) . Some recent works on mag- 
netic screening by a cc retion mounds (iPavne & Melatos 
20041: IPavne & Melatos] 120071 : I Vigelius & Melatod 12008 : 
Vigelius & Melatos! l2009h report that large mounds of mass 



~ 10 ~ Mq may form on the neutron star, which can then bury the 
field as the matter spreads on the surface. Ho wever several ques- 
tions regarding the ef fects of MHD instabilities dLitwin et al.l200ll : 
ICumming et al. 2001) remain to be addressed fully. Magneto- static 
solutions of accretion mounds have earlier b een found by sev- 
eral authors including Hameurv et al. dl983l). iBrown & Bildstenl 
dl998[) . IPavne & Melatos! d2004l) and iMukheriee & Bhattac harval 
d2012l) . It was shown in IMukheriee & Bh attacharvaT d2012h Qiere- 
after MB 12) that magneto- static solutions cannot be found for 
mounds beyond a threshold height (and mass), which may be 
indicative of the pres ence of MHD instabiliti es. Similar results 
were also reported in IPavne & Melatosl (120041) (hereafter PM04) 
where closed magnetic loops were seen to form beyond a threshold 
mound mass. 

In this paper we attempt to study the stability of the accretion 
mound by 2D axisymmetric M HD simulations with the PLUTO 
MHD code dMignone et alJ2007l) . The study of the full set of MHD 
instabilities in such mounds requires global 3D simulations. How- 
ever, results from 2D simulations would help to identify modes that 
grow despite of the restrictive assumption of axisymmetry. This 
will be a stepping stone to future 3D simulations where many other 
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modes may grow simultaneously. Here we investigate the pres- 
ence of interchange instabilities as predicted for such mounds by 
iLitwin et"aD ( 1200 lh . and also the physical cause of the threshold 
in mound mass obtained in MB 12. To study the latter, we add a 
small amount of mass to an existing GS solution and dynamically 
evolve the system to see if it settles to a new equilibrium state. This 
is carried out for different mound sizes up to the threshold mass, at 
which one expects MHD instabilities to be triggered if the threshold 
happens to be due to a physical effect. 

Our approach differs from that of PM04 in various aspects. 
We consider a cylindrical geometry with strict containment of the 
accreted matter in the polar cap, while PM04 consider spherical ge- 
ometry with mass loading on all field lines up to the equator. Also, 
we consider degenerate non-relativistic Fermi plasma near the po- 
lar cap surface instead of the isothermal equation of state used by 
PM04. As we consider densities as high as ~ 10 8 gcm -3 inside 
the mound, a degenerate non-relativistic plasma is more appropri- 
ate (see MB 12 for a discussion). 

Early models of accretion column formed by disc- 
magnetosphere interaction proposed hollow ring-like accre - 
tion column on neutr on s tar poles ( Basko & Sunyaev| dl976h . 
iGhosh & Lainbl dl978h and iGhosh & Lainbl dl979h ). Several au- 
thors have used hollow ring-like accretion col u mns t o fit the 
pulse p rofiles of HMXBs ( e.g. IShakura et all (fl99lh. iLeahvl 
dl99lh. iRiffert et all dl993h ). IPanchenko & Postnovl dl994l) and 



iKlochkov et al J <2008l) discuss effects of emission from two dis 
connected rings to explain shape of observed pulse profiles and 
nature of cyclotron features in the emission from Her X-l. Fol- 
lo wing the formalism of pulse profile decomposition developed 
by iKraus et al] dl995h . rins-like columns have been inferred for 
sources like Her X-l dKrausll200ll) 4U 1909+07 (lFurst et al.l201ll) 
A053 5+262 dCaballero etalfcOllh and V 0332+53 dFerrigno et al.1 
1201 lb . Even for LMXB sources, ring li ke polar cap models are pre- 
ferred for fitting pulse profiles (iPoutanen et al.ll2009t iKaiava etaD 
120111) . We therefore perform a study of the structure and stabil- 
ity of hollow accretion mounds to compare with results from filled 
mounds. We also perform simulations of CRSF emission from hol- 
low mounds, following the method described in MB 12. 

We structure the paper as follows: in Sec.|2]we outline the nu- 
merical set up involved in the problem. We discuss the solution of 
the Grad-Shafranov equation to determine the structure of the static 
mound. We also discuss details of the set up of the MHD simula- 
tions with PLUTO. In Sec. [3] we discuss the testing of the equilib- 
rium solution with PLUTO. In Sec. [4] we discuss the method and 
results of the perturbation analysis with PLUTO to investigate the 
stability of the mounds. In Sec. [5] we discuss the results of the sim- 
ulations of hollow mounds and we summarise the results in Sec. [6] 



2 NUMERICAL SET UP 

To test the hydromagnetic stability of the confined mound we first 
evaluate the equilibrium solution to the Magneto Hydrostatic equa- 
tions by solving the Grad-Shafranov (hereafter GS) equation. The 
solution of the GS equation is used as initial condition in PLUTO, 
where perturbation analysis is performed. In the following section 
we outline the solution of the GS equation and the set up of the 
simulation using PLUTO. 



2.1 Equilibrium solution from Grad-Shafranov equation 

For an axisymmetric system, one may write the magnetic field in 
terms of the flux function in cylindrical coordinates as 



VipxO 



(Be = 0) 



(1) 



Using eq. (Q} in the static Euler equation and using separation of 
variables in cylindrical coordinates using method of characteristics 
(as in MB 12) we get the GS equation for an adiabatic gas (p = 



A 2 ^ 



= -99- 



dZ 
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where g is acceleration due to gravity and density is given by the 
equation 
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Zo(V0 * s me mound height function which determines the shape of 
the mound. For our work we use the equation of state for a degen- 
erate non-relativistic zero temperature Fermi plasma with \i e — 2: 

^2 / n \5/3 



3.122 x 10^ 



fi e m p 
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(4) 



dynes cm 



10 6 g cm- 
Most of the mound will be dominated by degeneracy pressure ex- 
cept for a thin layer at the top (~ 4cm at IkeV plasma, see MB 12 
for a discussion). Thus effects of thermal stratification would play a 
limited role, and the zero temperature degenerate equation of state 
would be an adequate assumption. We solve the GS equation for 
an accretion mound of radius R p = 1km, on the poles of a slowly 
spinning neutron star of mass 1.4M© and radius R — 10km. The 
intrinsic field is assumed to be dipolar, which in the polar cap re- 
gion can be approximated as an uniform field along z (B p = B z). 
We consider Newtonian gravity with constant acceleration: 

-2 



; = -l 



x 10 1 



IAMqJ VlOkm J 
up is sim ilar to that in IHameurv et al.1 dl983h 



(5) 



Our set up is sim ilar to that in IHameurv et al.1 dl983|) and 
Litwin et al 1 booih . Following MB 12, we carry out most of our 
analysis for the mound height profile: 

where Z c is the central height of the mound and ip p = (l/2)BoR p . 
This is a smoothly varying parabolic profile in ip which describes 
a filled axisymmetric mound. We also discuss the GS solution for 
a hollow mound in Sec. [5] which is specified by the mound height 
function: 



Z c 

0.25 



0.25 



0.5 



(7) 



The GS is a coupled non-linear elliptic partial differential equation. 
We have solved the GS equation by an iterative under-relaxation al- 
gorithm with an inn er Successive Ov er-relaxation loop with Cheby- 
shev acceleration (iPress et al.l 1993) as is outlined in MB 12. For 
a given polar magnetic field (B p ), the solutions to the GS equa- 
tions are obtained up to a threshold height Z max , beyond which 
the numerical scheme does not converge to give an unique solu- 
tion. Details of the numerical algorithm and convergence of the GS 
solutions have already been discussed in MB 12. 
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z c 


B p 


N r x N z 


Al (Az ~ Ar) 


72m 


10 12 G 


1024 x 144 


~ 0.43m 


65m 


10 12 G 


1088 x 104 


~ 0.46m 


55m 


10 12 G 


1272 x 88 


~ 0.39m 


50m 


10 12 G 


1024 x 80 


~ 0.43m 


25m 


10 n G 


1920 x 72 


~ 0.2m 



Table 1. Sample resolutions for simulation runs. 



2.2 PLUTO setup: Initialisation 

We use the Godun ov scheme based MHD code PLUTO 
(Mig none et alj|2007b to test the stability of the confined mound. 
The solutions of the GS equation are used as initial condition in 
PLUTO. The GS solutions are imported into PLUTO using bi- 
linear interpolation. We use the MHD module of PLUTO to solve 
the full set of ideal magneto hydrodynamic equations: 





— + v • Vp + pV • v 

at 





(8) 


dv 


Vv + -B x (V x B) + - Vp 
P P 
OB 




(9) 




— + V x (v x B) 


- 


(10) 




— + v • Vp + pc s V • v 





(11) 



where the factor l/V^r is absorbed in the definition of mag- 
netic field and c 2 s is the speed of sound (which for adiabatic gas 
is c 2 s — jp/p). The system is closed by an equation of state (here- 
after EOS) which we choose to be either adiabatic (pe = p/(j— 1)) 
or barotropic for which p = p(p). In the second case eq. dTTb is re- 
dundant. To investigate the effects of pressure driven interchange 
modes and gravity driven modes, we perform perturbation analysis 
with the adiabatic EOS (see Sec.EQand Sec.0~2]>. PLUTO initial- 
isation and boundary conditions are provided in terms of primitive 
variables (p, v,p, B) defined in eq. (HJ) — eq. (UTT) . The computa- 
tion is carried out in conservative variables (p, pv, E, B), where 
E = pe + pv 2 /2 + B 2 /2 is the total energy density. 

We use the extended generalised Lagrangian multiplier 
(EGLM) scheme dMignone & Tzeferacosl J2010h , iMignone et al.1 
(12010 )) to preserve the V ■ B — constraint. The EGLM scheme 
preserves the divergence criterion by modifying the in duction equa- 
tion (eq. [Tot with a scalar field function ^glm foedner et al.l 
120021) and also the energy momentum equations with extra source 
terms. This scheme transports the non-zero divergence errors to the 
boundary of the domain at the fastest possible characteristic speed, 
and damp them at the same time. 

For our problem, we have foun d that the HLL Riemann so lver 
dTorol2008h . HLLD Riemann solver jMivoshi & Kusancl2005h and 
TVD Lax-Friedrichs solver dTorol l2008h combined with EGLM 
scheme provide solutions free from numerical instabilities. Due 
to the presence of very sharp gradients in the physical quantities, 
higher order schemes need to be employed to reduce numerical 
errors. A third order Runge-Kutta scheme is used for time evolu- 
tion and a third order accu rate piece- wise parabolic int erpolation 
scheme (PPM scheme as in Colella & Woodward (1984)) has been 
employed. 

The simulations were set up using square cells (Ar ~ Az) 
to minimise numerical errors. The resolutions used were less than 
~ 0.5m as listed in Table. Q] for some sample runs. The physical 
variables in PLUTO are scaled to non-dimensional forms before 



Field line plot for Z c = 65m 




Radius in Km 

Figure 1. Field lines for a mound of height Z c = 65m with polar un- 
loaded field Bp = 10 12 G. The dash-dotted line in red denotes the top of 
the mound beyond which density is zero. The total mass of the mound is 
~ 1.63 x 10 _12 M©. The dashed blue box in the middle is the PLUTO 
computation domain, chosen to keep Alfven velocities non-relativistic. The 
range of density is ~ 2.1 x 10 6 - 6.7 x 10 6 g cm -2 at the top of the 
mound and ~ 3.02 x 10 7 — 5.7 x 10 7 g cm -2 , at the bottom. 



initialisation. For example for mounds with polar magnetic field 
Bp = 10 12 G, we use p = 10 6 gcm -3 as the density unit, L = 
10 5 cm as the length unit, Bo — 10 12 G as the magnetic field unit 
and Vao = B I^J^p = 2.82 x 10 8 cm s _1 as the velocity unit. In 
these units, time is measured in units of £a = Lo/Vao = 3.55 x 
10 _4 s, which can be taken as the mean Alfven time, while the scale 
velocity is the mean Alfven velocity. An unique Alfven velocity 
cannot be prescribed for the whole domain as the Alfven speeds 
will vary over the domain depending on local density and magnetic 
field. 



2.3 Boundary Conditions 

For stability studies, we run the simulations with either fixed 
boundaries where quantities are kept fixed to initial values (Q = 
Qo) or fixed gradients where the initial gradients are preserved. 
The fixed gradient boundary implies outflow of perturbed quanti- 
ties as gradients of perturbations are set to zero (VQ = VQo + 
VQ — ► VQo,VQ = 0). The standard outflow boundary condition 
(VQ = 0) is inapplicable for our problem as the initial solution 
has non-zero gradients at the boundaries of the domain. The fixed 
gradient boundary condition is applied to the upper and the right- 
most boundary. For filled mounds, the inner left boundary is kept 
fixed as it is close to or equal to the axis of the column. For hollow 
mounds, the inner left boundary is kept at a fixed gradient to allow 
for inward flow of perturbed matter. The bottom boundary is kept 
fixed to simulate a hard crust. The set-up with fixed gradients on 
the outer sides and fixed crust gives numerically stable solutions, 
as tested from simulations of the equilibrium solutions obtained 
from GS-solver (see Sec.|3}- 



3 EQUILIBRIUM STUDIES 

The GS solutions for adiabatic mounds have density profiles which 
go to zero beyond Zo(ip) (see eq. $3$). To avoid unrealistic Alfven 
velocities, we restrict the computation domain inside the mound 
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Figure 2. Energy components for zero-mean random perturbation run, nor- 
malised to their initial value. Magnetic energy is normalised to 3.7 x 10 22 , 
internal energy to 8.9 x 10 23 ergs and gravitational potential energy to 
6.7 x 10 23 ergs. The internal and gravitational energy components remain 
almost constant (~ 0.02% change from initial value). The magnetic energy 
initially decreases as the pockets of perturbed matter settle down, eventu- 
ally returning to its initial value. This indicates that the system is stable, and 
when perturbed, settles to a energy state close to the original equilibrium 
value. 



such that Alfven speeds in the mound are non-relativistic. A typi- 
cal computation domain is depicted in Fig.Q]for a mound of height 
Z c — 65m. We first evolve the initial equilibrium solution without 
applying perturbation in order to check the stability of the numer- 
ical schemes and also to study the effects of initial transients con- 
tributed by the numerical errors accumulated in interpolating the 
solution from GS grid to PLUTO domain. 

The solutions have been evolved to £ ^ 80£a for different 
choices of schemes. For the set of schemes outlined in Sec. l2.2l and 
Sec. 12.31 the equilibrium solution remains intact, with very small 
build up of internal flow velocities. For example, for a mound of 
height Z c — 72m, at t ~ 80£a, the maximum velocity is ~ 7.5 x 
10 -4 in normalised units (~ 2.15 x 10 5 cms _1 , which is much 
smaller than typical scale velocities). This shows that the schemes 
used are free from artificial numerical effects and also verifies the 
validity of the equilibrium solution obtained from the GS solver. 



4 PERTURBATION ANALYSIS 

We perturb the equilibrium solution by adding a normalised pertur- 
bation field £(r, z) to any of the physical quantities 

Q = Qo(l + rt(r,z)) (12) 

where r\ is a positive number signifying the perturbation strength. 
The perturbations are kept away from the boundaries on all sides. 
This is to preserve the equilibrium at the boundary layers and avoid 
spurious interaction with the boundary. For our studies we apply a 
random perturbation on the density inside the simulation domain, 
namely £ is assigned a random value at each grid point within the 
perturbation zone. The edges of the perturbing region are smoothed 
with an exponential function to avoid sharp gradients which can 
lead to spurious effects. The lack of any preferred perturbation scale 
should allow the growth of the fastest growing modes. The pertur- 
bation analysis is performed for mounds of different heights up to 



Figure 3. Energy components for random positive perturbation run (77 = 
5%), normalised to their initial value. Magnetic energy is normalised to 
5.2 x 10 22 ergs, internal energy to 9.9 x 10 23 ergs and gravitational po- 
tential energy to 7.5 x 10 23 ergs. The initial energy is dominated by internal 
and gravitational energy. The gravitational and internal energy decrease as 
the system move to a lower energy state following the perturbation. The 
magnetic energy is seen to increase due to stretching of field lines due to 
internal flows. 

the threshold height Z max beyond which the GS- solver does not 
converge, as has been found in MB 12. 

4.1 Zero-mean perturbations: interchange modes 

Zero-mean random perturbation with (£) = 0, implies rearrang- 
ing of density from the equilibrium solution without adding any 
net mass. In this case, the system quickly converges to stable 
pockets of perturbations, irrespective of perturbation strength (77 in 
eq.[l2]). See Fig.|4]for the results of a run with perturbation strength 
77 = 10%. The system settles down to an energy state close to the 
original equilibrium value (see Fig.|2j. However, for larger pertur- 
bation strengths, a longer time is taken to relax into stable pockets 
of perturbed matter. For example a mound with B p = 10 12 G and 
Z c = 65m stabilizes after t ~ 1£a for 77 = 2% and t ~ At a for 
77 = 10%. 

The perturbation tests have been carried out for mounds of 
different heights and polar magnetic field strengths. No instabili- 
ties are seen at the threshold mound heights, e.g. Z c ~ 72m for 
B = 10 12 G and Z c ~ 25m for B = 10 n G etc. The simulations 
show that the mounds are stable with respect to small departures 
from equilibrium resulting from rearrangement of flux tubes . Thus 
interchange or ballooning modes are not seen in 2-D axisymmetric 
simulations of the mounds. 



4.2 Adding excess mass to equilibrium solution 

In order to study the effect on the mound of the addition of matter 
which eventually descends due to gravity, we apply a positive defi- 
nite random perturbation field: (£) > on the density without any 
corresponding change in pressure. Such a change in density implies 
local departure of A; a d from that in eq. |4] In this work we do not at- 
tempt to model the exact composition of the accretion mound. The 
perturbations were set up to ensure that the added matter is heavier 
than its surroundings and will descend due to gravity, thus trigger- 
ing the gravity driven modes. However, a change in k a d can indeed 
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0.4 0.5 0.6 0.7 



Figure 4. Over-density: (p — p e q ) /peq for zero-mean perturbation runs for a mound of height Z c = 65m, polar magnetic field B p = 10 12 G and perturbation 
strength 77 = 10%. p e q is the unperturbed density from the equilibrium solution. The vertical axis is the height above neutron star surface in kilometres. The 
horizontal axis is the radius (cylindrical geometry) in kilometres. The PLUTO simulation was carried out with a grid of size 1024 x 120. Random perturbation 
is provided within a rectangular box inside the domain, away from the boundaries. The edges of the perturbation region are smoothed exponentially. The 
perturbation slowly weakens and relaxes into stable pockets of perturbed density by t ~ At a (bottom panel). The magnetic field lines are plotted in black. 




Over density at t = 50.0000 t A 




Figure 5. Over-density: (p — p e q)/peq at different times for a positive density perturbation with strength 77 = 3% in a mound of height Z c = 65m and 
polar magnetic field B p ~ 10 12 G. The simulation was carried out with a grid of size 1088 x 104. Horizontal and vertical axes are the same as in Fig. [4] The 
perturbations result in the formation of closed loops but the solution eventually settles down to a steady state. 




B/B eq at t = 37.0000 t A 




0.4 0.5 0.6 0.7 



Figure 6. Magnetic field magnitude normalised to the local equilibrium value for the simulation described in Fig. [5] Bunching of field lines forms pockets of 
excess field over equilibrium value, which eventually get smeared and start to dissipate. 
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Over density at t = 18.2500 t A 




Over density at t = 38.5000 t A 




Over density at t = 50.0000 t A 




Figure 7. Over-density (p — p e q)/peq at different times for a positive density perturbation with strength 77 = 5% in a mound of height Z c = 65m and 
Bp = 10 12 G. The simulation was carried out with a grid of size 1088 x 104. Horizontal and vertical axes are the same as in Fig. [4] Reconnection of field 
lines forms closed loops at multiple sites. The system does not relax to any steady state solution within the duration of the run. The closed loops grow with 
time indicating the onset of unstable modes. 



B/B eq at t = 18.2500 t A 




0.4 0.5 0.6 0.7 



B/B eq at t = 38.5000 t A 

I ^ ^ ^ ^ 1 H3.73 

°- 04 r ^ ^ 2.80 
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B/B eq at t = 50.0000 t A 




0.4 0.5 0.6 0.7 



Figure 8. Magnetic field magnitude normalised to the local equilibrium value for the simulation described in Fig. [7] Bunching of field lines cause pockets of 
excess field over equilibrium value which do not settle to any steady state. 

occur due to changes in chemical composition e.g. 77 ~ 5% local 
perturbation on a Z c ~ 65m mound would correspond to a change 
of mean molecular weight by A/i e ~ 0.1. 

The added mass settles down along the field lines, dragging 
and distorting the equilibrium field configuration in the process. 
For small perturbation strengths (77 ~ 1% for mound of height 
Z c = 65m) the matter quickly settles down to a new equilib- 
rium, without appreciable distortion of the field lines. With an in- 
crease in 77 beyond a threshold, e.g. t\t ~ 3% for Z c — 65m and 
Bp — 10 12 G mound, magnetic Rayleigh-Taylor type instabilities 
are triggered by the descending heavier matter and results in the 
formation of closed loops due to reconnection of field lines (see 
Fig. 13 Bunching of field takes place in the radial direction (e.g. 
Fig-EJ and the system eventually relaxes to a steady state. 

Further increase in perturbation strength, e.g. 77 ~ 5% for 
Z c — 65m, disrupts the equilibria completely. Several closed loops 
are formed across the perturbed region (see Fig [7] and Fig. |§J. In- 
dividual closed loops merge to form larger knots without showing 



Note that although the simulation is ideal MHD, numerical resistivity 
allows dissipation and reconnection to occur. 



any signs of decay. From Fig. [3] we see that gravitational potential 
energy and internal energy decreases from initial value, whereas 
magnetic energy increases with time. This indicates that internal 
flows stretch and twist the field lines converting internal energy and 
gravitational energy to magnetic energy. The system does not relax 
to a steady state within the run time of the simulation (t ~ 50£a). 
Thus for a mound with Z c = 65m and B p — 10 12 G, the thresh- 
old perturbation strength is t\t ~ 3% beyond which gravity and 
pressure driven modes disrupt the MHD equilibria. 

Convergence has been tested by running the simulations for 
successive higher resolutions: .e.g. for Z c — 65m, B p = 10 12 G 
with positive random perturbation of strength 77 = 5% simulations 
were carried out for resolutions (1088 x 104), (2176 x 208) and 
(4352 x 416). It was seen that MHD instabilities persist on in- 
crease of resolution. Increase in resolution reduces numerical re- 
sistivity, thus decreasing cross field diffusion. The field lines are 
then more prone to be deformed by gravity driven modes triggered 
by the weight of the overlying matter. 

With an increase in mound height, it is easier to excite such 
unstable behaviour. The threshold perturbation strength is larger 
for mounds of smaller height: for Z c = 45m and B p = 10 12 G, 
rjT ~ 7%. Mounds near the GS threshold height Z max (~ 72m for 
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Hollow mound of height Z c = 45m 




0.0 0.2 0.4 0.6 0.8 1.0 
Radius in Km 



Figure 9. The field lines from GS solution for a hollow mound with mound 
height function given by eq. {7}, Z c = 45m and B p = 10 12 G. The max- 
imum height Z c occurs at ~ 698m from the axis. The red-dashed line 
represents the top of the mound. 



CRSF for parabolic mound Z c ~ 45m CRSF convolved with Gaussian 
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Figure 10. The magnetic field at the top of the hollow mound in Fig. [9] 
Field lines are pushed on either side of the apex (r ~ 698m) of the mound 
resulting in decrease in field at the apex and increase in field strength on 
either side. Starting from a polar magnetic field strength B p = 10 12 G, 
from our GS solution we get minimum field at the top ~ 6.63 x 10 n G and 
maximum field of ~ 2.33 x 10 12 G. 



B p = 10 12 G; - 25m for B p = 10 n G) are only marginally stable 
at t\t — 1%. Thus, mounds higher than a threshold (as previously 
obtained in MB 12) are prone to gravity driven Rayleigh-Taylor and 
pressure driven instabilities on addition of excess mass, and stable 
magneto- static solutions cannot be obtained. 

5 HOLLOW MOUND 

5.1 Grad-Shafranov for hollow mounds 

For systems with magnetospheric accretion, mass loading at the 
accretion disc takes place over a finite rang e of accretion disk radii 
(A r ~ 0.03#a, Ra = Alfven radius e.g. [Ghosh & Lambl (il978l. 
Il979l) ). The inner edge of the polar cap ring □ for such systems will 



2 which corresponds to the outermost radius in the accretion disc ~ Ra + 
A r , where mass loading begins. 



Figure 11. Top: CRSF from a filled mound of central height Z c ~ 45m 
and Bp ~ 10 12 G. The right panel gives the spectra convolved with a 
Gaussian (standard deviation 10% of local energy) to simulate finite de- 
tector resolution. Bottom: CRSF from a hollow mound with Z c = 45m 
and Bp = 10 12 G, with the right panel giving the convolved spectra as be- 
fore. The CRSF from hollow mounds show a much broader spectra due to 
contribution from different parts of the mound with large variations in the 
magnetic field. 



be 

^ = ^( 1 "^) (13) 

while the outer edge of the polar cap radius is (R s / Ra) 1 ^ 2 R s 
jPoutanen et alj|2009k R s being the neutron star radius. For small 
values of A r the columns would be hollow and thin walled. On 
the surface of the star this would create an accretion ring around 
the polar cap instead of a filled mound. To model such an accretion 
ring, we choose the mound height function to give a hollow mound 
in which the density falls off to zero both at the axis and at the polar 
cap radius. 

For the solution presented in Fig. [9] we use a mound height 
profile as in eq. (0 with Z c — 45m and B p = 10 12 G. The solution 
shows considerable distortion of field lines on both sides of the 
apex (r ~ 698m). This is in contrast to the case of filled mounds, 
where curvature of field lines occur towards the outer edge. Larger 
curvature of field lines allow larger mass to be accumulated per 
flux tube, as compared to that of filled mounds. Hence, although 
the central part is hollow, the total mass contained in the hollow 
mound (M ~ 5.87 x 1O _13 M ), is comparable to that of a filled 
mound of the same height and field ( M ~ 5.09 x 10~ 13 M© for 
Z c rsj 45m and B p ~ 10 12 G and a parabolic profile as in eq.|6j. 

The family of GS solutions for hollow mounds behave simi- 
larly as for filled mounds. With increase in maximum mound height 
Z c , the GS solutions show larger curvature of field lines on both 
sides of ridge apex. The GS solutions fail to converge for mounds 
greater than a threshold height for a given magnetic field. For the 
mound height profile of eq. (0, the threshold height is around 
Zmax ~ 47m for a polar magnetic field B p = 10 12 G. 
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Figure 12. Over density (p — p e q)/peq and field line for hollow mound of maximum height Z c = 45m and polar magnetic field B p = 10 12 G, with a 
positive density perturbation of strength rj = 5%. The simulation was carried out for a grid of size 1144 x 136. The vertical and horizontal axes are the same 
as in Fig. |31 The perturbation results in formation of closed loops at multiple sites near the centre, very early in the simulation run. 



B/B pn at t = 1.80000 t A 



0.02 



0.01 



0.02 







0.60 0.65 

B/B eq at 


0.70 0.75 0.80 
t = 9.15000 t A 


[ 4 < 


• ► * ) 


0.60 0.65 0.70 0.75 0.80 

B/B eq at t = 30.0000 t A 







0.60 



0.65 



0.70 



0.75 



0.80 



2.22 
1.67 
1.11 
0.56 
1 0.00 



2.53 
1.90 
1.27 

■ 0.63 

■ o.oo 



1.89 
1.44 
0.99 

-J 10.53 

0.08 



Figure 13. Magnetic field magnitude normalised to the local equilibrium value for the simulation described in Fig |12l Bunching of field lines in radial direction 
causes alternate regions of enhanced field strengths. The closed loops and pockets of enhanced fields migrate to the radial boundaries and eventually dissipate. 



5.2 Stability analysis of Hollow mounds 



Using the GS solutions for hollow mound, we perform stability 
analysis with PLUTO. The results are similar to that of a filled 
mound. Zero-mean density perturbations do not show growth of 
the perturbed region, indicating that the mounds are intrinsically 
stable with respect to interchange modes. For positive perturba- 
tions in density, closed loops are formed after a threshold pertur- 
bation strength. See Fig. [12] and Fig. [13] for the results of a run 
with r\ — 5%. The closed loops form quickly within a few Alfven 
times and migrate away from the center, on both sides of the cen- 
tral height. This results in the formation of alternate regions of en- 
hanced and reduced magnetic field due to the bunching of field 
lines, which have considerable departure from equilibrium solu- 
tion. The field knots dissipate gradually as they migrate outwards. 



5.3 Cyclotron lines from hollow mounds 

Following the algorithm outlined in MB 12, we have simulated the 
cyclotron resonance scattering features (hereafter CRSF) that will 
be observed in the emitted spectrum from a hollow mound. The 
spectra have been calculated by integrating the emission from dif- 
ferent parts of the mound towards a given line of sight (hereafter 
los). We assume a Gaussian absorption profile whose depth and 
width are evaluated from interpolated results of ISchonherr et al.l 
d2007h for the slab 1-0 geometry. As in MB 12, the line centre of 
the CRSF is obtained from the expression 

E n = nEcoVT^l (l-^ ( 51 f^ v ) sin 2 6u) (14) 

where n = 1, 2, 3... is the order of the harmonic, E c o = II.6B12 
in keV, O^b is the angle between the direction of emission and local 
magnetic field and u = r s /r, r s being the Schwarzschild radius. 
Emission from the inner part of the hollow mound may be blocked 
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Plasma (3 for a mound of height Z c ~ 65m 




Figure 14. Plasma (3 (ratio of plasma pressure to magnetic pressure) for a 
GS solution of a mound of height Z c ~ 65m and B p ~ 10 12 G. The verti- 
cal and horizontal axes are the height and radius respectively, expressed in 
kilometres. The maximum plasma (3 (~ 911) occurs along the central red 
horizontal patch near the regions of maximum curvature of the magnetic 
field lines (represented in white). At the regions of high (3, the plasma is pri- 
marily supported by tension from curvature of field lines. Such regions are 
prone to pressure driven instabilities, and show formation of closed loops 
when perturbed. 

by the walls on the opposite side. In Appendix. (A]) we explain the 
scheme we follow to account for such shielding. 

For the simulated spectra shown in Fig. Qj] we consider emis- 
sion from a single pole at inclination angle r\ v — 10° and a los 
at i = 60°, both measured from the spin axis. The spectrum 
shows multiple absorption features due to the large variation of field 
strength at the top of the mound (see Fig.[l0j. The different absorp- 
tion features correspond to emission from different locations on 
top of the mound, with different magnetic field values. The nature 
of this spectrum is significantly different from that expected from a 
filled parabolic mound of the same height (see Fig.[TT). When con- 
volved with a Gaussian of standard deviation ~ 10% of the local 
energy, to simulate the finite resolution of a detector (see MB 12 for 
details), the spectrum becomes a broad absorption feature. 



6 DISCUSSION AND SUMMARY 

(i) Absence of interchange mode instabilities: In this paper 
we have tested for the stability of magneto- static accretion mounds 
by MHD simulations using the PLUTO MHD code. From perturba- 
tion analysis we conclude that mounds are stable with respect to in- 
terchange or ballooning mo des in 2-D ax i symm etric simulation^. 
Linear stability analysis by iLitwin et"al . (2001) predict the onset 
of ballooning modes for a threshold plasma (3 ((3 = p/(B 2 /8tt)). 
However such modes are inherently multi-dimensional in nature, 
with finite toroidal and z ero poloidal wave vectors, normal to the 
local magnetic field (see Freidberd dl982h for a review of MHD 
instabilities in confined plasma). Hence such modes cannot be ex- 
cited in an axisymmetric 2D simulation. 

Litwin's approximate analytical estimates give a threshold /3t ~ 
11.7(Rp/Z c ) for 7 ~ 5/3, beyond which MHD instabilities will 
set in. For GS solution of a filled mound with Z c = 45m and 
Bp = 10 12 G, we get maximum /3 ~ 293, which is close to 

3 Note that in this paper we consider a T = 0°K Fermi gas. How- 
ever, finite plasma t emperature can induce additional thermal modes 
( Cumming et al. 2001) which have not been explored here. 



Litwin's threshold for the same mound /3t ~ 260. For higher 
mounds, /3t decreases with increase in Z c , and is much smaller 
than the maximum (3 obtained from our GS solutions. For example, 
for a filled mound with Z c = 65m and B p = 10 12 G, f3 T ~ 180 
whereas maximum f3 ~ 911 from the GS solution (see Fig IT4t. 
Hence results from 2-D simulations cannot rule out the presence 
of su ch modes i n a 3-D set up. Also, interchange mode instabili- 
ties faienll 19841) can be excited in 3D simulation runs, as is seen in 
other examples of confined plasmas, e.g. in tokamak reactors. Work 
on 3-D stability analysis of accretion mounds is currently underway 
and will be addressed in a forthcoming publication (Mukherjee, 
Bhattacharya and Mignone in preparation). 

(ii) Instabilities due to excess mass: From our 2-D simula- 
tions we have found that addition of excess mass destabilizes the 
equilibrium due to gravity driven magnetic Rayleigh-Taylor type 
instabilities. For mounds with higher mass, the GS solutions have 
large radial (horizontal) component of magnetic field, which be- 
ing p erpendicular to gravity are also prone to Parker ty pe instabil- 
ities JCumming et all 1200 ll : iMelatos & Phinnevl 120011) . Topologi- 
cally disconnected closed loops are formed beyond a threshold per- 
turbation strength tjt. 

Fro m the expression o f the energy integral for linear perturba- 
tions dLitwin et alj|200ll) on an adiabatic plasma (p = kp 1 ), we 
have 

5W=y d\{^ + g (V • £ ± + 2n c ■ £ ± ) 2 

+ 7P(V.£-2K 9 .£) 2 tl5) 
- 2(k c + V<j>/(2c 2 s )) ■ $± (Vp + P V<P) ■ £ ± } 

where £ is the plasma displacement, B = Vx x B) is the per- 
turbed magnetic field, n c = (b- V)b is the magnetic field curvature 
vector, c s is the sound speed and (j) the gravitational potential. B^ 
is zero for our case. 

Instabilities will develop if the negative contribution from any 
(or all) of the terms containing field curvature, pressure gradient 
and gravity overcomes the stabilizing effects of the magnetic and 
pressure compression terms. Q Hence, it is not a surprise that the 
closed loops are formed in regions with the largest curvature in field 
lines. This also corresponds to the regions with high plasma ft, e.g. 
the red region in the middle of Fig. [14] where f3 ~ 911. Pressure 
driven instabilities typically lead to a threshold plasma (3 beyond 
which instabilities are triggered (e.g. Freidbergl (ll982klLitwin et aD 
For mounds near the stability threshold, e.g. Z c ~ 72m at 
B p — 10 12 G, the maximum plasma (3 is as high as ^ 1.26 x 10 4 . 

The magnitude of t\t decreases with increase in mound height, 
with rjT — > as Z c —> Z max , indicating inherent unstable nature 
of the mound for the modes under investigation. This corroborates 
the result of MB 12 that GS solutions do not converge beyond a 
threshold height. The tests involving addition of mass are not meant 
to reflect realistic accretion rates. Although the amount of excess 
mass added in our simulations is small (~ 7.6 x 10 _15 M© for 
rj = 5% perturbation on 65m mound), in a real system such mass 
will be accumulated slowly as mounds of larger mass are built. Ef- 
fects of such inflow of material on an initially static mound have 
not been addressed here. However, from our current 2-D simula- 
tions we conclude that for large mound masses, gravity and pres- 
sure driven modes result in the onset of MHD instabilities and no 
static equilibrium solution can be found beyond a threshold Z max . 

4 Necessary and su fficient condition for instability is: 5w < 
teernstein etaO 19581) . 
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Magnetic field and velocity unit vectors at t ~ 5t A 




Radius in km 



Figure 15. Ratio of the magnetic field to its local equilibrium value for a barotropic simulation with random velocity perturbation of strength r\ = 15% of local 
sound speed (maximum initial velocity ~ 0.35, in normalised units). The velocity unit vectors are plotted to show the nature of the flow. Bunching of magnetic 
field takes place in the radial direction as local eddies are formed. The system settles down to a steady state with flow velocities less than ~ 7.54 x 10 -3 (in 
normalised units) at t ~ 5t a ■ 



Buoyancy related instabilities due to the formation of topolog- 
ically disconnected closed loops hav e previously been rep orted 
in the static mound si mulations of iHameurv et alj |l983) and 
Payne & Melatos ( 2004), and also dynamic MHD simulations by 
Vigelius & Melatosl (2008) (hereafter VM08). However the thresh- 
old mass of the mound for the formation of closed loops in PM04 
and VM08 is M ^ 10 _5 M©, which is much larger than the mass 
of the mounds in our present work. This may be due to the follow- 
ing differences in approach: 

(a) PM04 and VM08 in their treatment consider spherical 
polar geometry and populate all field lines up to the equator, 
whereas we confine the accretion mounds strictly within the po- 
lar cap radius. Populating all field lines up to the equator pro- 
vides lateral pressure support to the polar mound which can then 
hold a larger mass. 

(b) Plasma pressure due to isothermal EOS by PM04 and 
VM08 is several orders of magnitude less than the degenerate 
Fermi pressure used in our treatment, which results in higher 
plasma f3 in our simulation. Such a system is more p rone to pres- 
sure driven MHD instabilities e.g. Freidberg(1982). 

(iii) Adiabatic vs barotropic: We have also performed 
barotropic simulations with PLUTO for which the energy equation 
becomes redundant as pressure is evaluated from p = kp 1 , with 
k a constant. This is similar to the isothermal set up of MHD sim- 
ulations. Results from adiabatic and barotropic modes are similar 
when perturbations are applied to velocity and magnetic fields. See 
Fig -Elf or the results of velocity perturbation with barotropic simu- 
lation (77 = 15% of local sound speed). The magnetic field bunches 
in the radial direction and local eddies are set up. The system settles 
down to a steady state with flow velocities reduced by more than 3 
orders of magnitude at t ~ 5£a- Similar results are also obtained 
for adiabatic EOS. 

However density perturbations behave differently in barotropic 
and adiabatic simulations. For a barotropic simulation, positive 
density perturbations on an initial static equilibrium create regions 
of excess pressure. The perturbed regions with high local pres- 
sure overcome the downward gravitational force and are quickly 
transported vertically upwards. Hence to study the effect of gravity 
driven modes due to the descent of added matter, adiabatic simula- 
tions have been performed in this work. 

(iv) Hollow mounds - structure and stability: We have solved 
the GS equation for mounds with hollow interiors. The hollow 
mounds show considerable distortion of the magnetic field on both 
sides of the maximum height to support the confined matter. There 
is a decrease in field near the ridge apex as field lines are pushed 
to either side. Closed loops form when excess mass is added to the 
equilibrium solution. The closed loops migrate to either side and 
eventually dissipate. 



The fixed gradient boundary condition can induce artificial sta- 
bility as it results in line tying type boundaries, which are known to 
give extra stability. In a real system, the plasmoids will be eventu- 
ally ejected from the system. Plasma travelling inwards on closed 
loops may then eventually fill up the hollow. However there was no 
significant mass loss seen in our 2-D simulations. 

(v) Hollow mounds - CRSF: CRSF from hollow mounds have 
been explored. From the simulation of the spectra integrated over 
the entire mound we see that: 

• Cyclotron emission from the top of hollow mounds show 
complex fundamental features in the line shape (harmonics have 
not been evaluated), due to the large variations in magnetic field 
on top of the such mounds. This is similar to what is observed in 
the sp ectra of V0332+53 (iMowlavi et alJbOOd iNakaiima et all 
2010) which is conje ctured to have a hollow column geometry 
dFerrigno et al .11201 lh . Complex line shapes have also been pre- 
dicted previously for strong non-dipolar local magnetic field by 
lNishimural(l2008Ll201lh . 

• Convolving the CRSF with a Gaussian to account for finite 
energy resolution of detectors, we see that the resultant CRSF 
has the structure of a broad absorption envelope. 

Thus CRSF from hollow mounds will be characterised by broad 
line widths and complex structures in the line shape, which may be 
observed with improved detector resolution. 

Thus we conclude from this work that accretion mounds on 
neutron stars in HMXB are stable up to a threshold height and mass, 
beyond which MHD instabilities will disrupt the equilibria. Struc- 
ture and stability of hollow mounds have been explored. It is shown 
that CRSF from such mounds will be characterised by broad fea- 
tures with a complex line shape. More work needs to be done to 
explore the 3-D stability of such systems and the effect of non- 
axisymmetric modes on the field structure and cyclotron emission 
from such mounds. 
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Figure Al. Top: A 3-D schematic representation of the hollow mound. The 
vectors and r c — r denote the plane where the path of the emitted ray 
to the observer lie. Bottom: a cross section of the mound along the plane 
of the emitted ray to the observer, and the location where the plane cuts the 
mound on the opposite side. 

ISSI, Berne and discussions with the Magnet collaboration which 
have benefited the paper. 



APPENDIX A: SHIELDING OF RADIATION FROM 
INNER WALLS OF HOLLOW MOUNDS 

In HMXBs an accretion column is formed by the infalling mat- 
ter after it passes through a shock which may be several kilo- 
metres from the surface of the star, d epending on the accre - 
tion rate (e.g. iBasko & Sunvaevl dl976h . IBecker & Wolffl (l2007h . 
iBecker et alT(l2012l) ). In this work we consider the spectra gener- 
ated from the mound without incorporating the effects of scatter- 
ing from the overlying accretion column. This is valid for systems 
with low accretion rates and optically thin columns. The emission 
from the mound will then be directly visible and effects of overly- 
ing column will be small. However for systems with optically thick 
columns and large accretion rates, the emission from the mound 
will be obscured by scattering and absorption in the column. A 
proper Monte-Carlo simulation of the radiative transfer through the 
column must be carried out to address such cases, which will be 
reported in a future work (Kumar, Bhattacharya and Mukherjee in 
preparation). 

The rays of light coming from the hollow region can be 
blocked by the inner walls of the mound on the opposite side. Such 
rays will not contribute to the total spectra. The path of the emitted 
ray lies in the plane defined by the radius vector from the origin 
(centre of the neutron star) to the point of emi ssion (r) and t he uni t 
vector along the line of si g ht (see f or e.g.lBeloborodovl d2002h. 
Pout anen & Beloborodovl d2006h and iMukheriee & Bhattacharval 
120121)). To exclude rays that may be blocked by the inner walls 
of the hollow mound, we first find the point where the plane de- 
fined by r and passes through the top of the mound r c as in 
Fig. lAll The radial and vertical coordinate of r c (r c and z c respec- 
tively) are found by fitting a polynomial to the top of the mound 
obtained from the GS solution and evaluating the coordinate where 
z is maximum. Since the three vectors r, and r c lie in the same 
plane, the angular coordinate (j) c of r c is found from the condition 

r c • (r x n^) = (Al) 



Following MB 12 we use the following definitions for the vec- 
tors = (n^ x ,n^ y ,n^ z ) = (sin i sin cj, sin i cos a;, cos i) and 
r = (x, y, z) = {p cos (f), p cos T] p sin^ + (£ + R s ) sin?7p, (£ + 
R s ) cos rj p — psin rj p sin 0} where i is the azimuthal angle of the 
observer's line of sight with respect to the spin axis, uj is the spin- 
phase angle, (p, 0, £) are coordinates of the emitting region in the 
polar cap frame with cylindrical coordinate system, R s is the neu- 
tron star radius and r\ v is the azimuthal angle of the centre of the 
polar cap. Using the above, we can rewrite eq. (IA11 as 

A c cos (f) c + B c sin C + C c = (A2) 

where 

A c = p c (yn^ z - zn^y) 

B c = p c cos Tj p (zn^x - xn^z) - p c sin r\ v (xn^ y - yn^ x ) 
C c = (£ c + R s ) {sin rj p (zn^ x — xn^ z ) 
+ cos Tjp(xn^ y - yn^x)} 

Eq. IA2I is solved using a modified Newton-Raphson scheme fol- 
lowing IPresTetal] |l993). After finding the coordinate of r c we 
evaluate the angle 0i c (see Fig. lAll) between the local normal (hi) 
and the radius vector from the point of emission to the top of the 
mound on the other side. 

cosft c = n r i rc ~ r , (A3) 
|r c -r| 

The normal vector is found as outlined in MB 12 by evaluating the 
slope m s = d^top/dp of the function ^ top = /(p) (p being the ra- 
dial coordinate) that fits the top profile of the mound obtained from 
the GS solutions: hi = { — sin S cos 0, — sin S cos r\ v sin + 
cos S sin ri P , cos S cos rj p + sin S sin rj p sin 0}, where sin S — 

, ms and cos S = , 1 = . Using the above definitions of the 
V 1+m s v 1+m s 
vectors, one can write 

nj • (r c — r) = cos# s (£ c - £) + sin # s (p cos - p c cos0 c ) 
x (sin <f> + cos <j)) 
|r c -r| 2 = p*+pl + (t-t c f 

— 2pp c (cos cj) cos (j> c + sin cj) sin <j> c ) 

Any ray with emission angle larger than 6i c will not contribute to 
the total spectra. This implicitly assumes that light will travel in a 
straight line and curvature effects from bending due to gravity is 
ignored for such short paths. 

More accurate methods should be used to calculate the tan- 
gent vector from the point of emission to the mound surface on the 
other side. However, this involves more computation, and for sharp 
profiles of the hollow mound used, approximating the tangent point 
as the top of the mound will result in only a small correction. 
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